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A PUNCHED CARD METHOD FOR SUCCESSIVE INTERVALS SCALING^ 



Samuel J. Messick and Ledyard R. Tucker 
Princeton University and Educational Testing Service 

and 



Harry W. Garrison 
Educational Testing Service 



INTRODUCTION 

The scaling method of successive intervals has appeared in several forms 
(e.g., 1, 3, 5) with variations In the computational labor required depending 
upon the simplifying assumptions made . Generally, the more rigorous the 
solution, the more laborious the computations (see reference 4); but, in any 
event, as the number of stimuli to be scaled increases, the amount of computation 

required by any suggested procedure soon becomes prohibitive . Recently, 

2 
however, an iterative least squares solution to successive intervals was 

developed, for which punched card procedures were appropriate; and thus a 

feasible computational routine for scaling large numbers of stimuli was made 

available without sacrificing the rigor of a least squares approach. 

A punched card procedure for a general weighted least squares solution 

to successive intervals is described below. Use of zero weights permits 

application of this procedure to cases of incomplete data . The procedure 

outlined suggests use of an IBM 101 Electronic Statistical Machine for obtaining 

frequency counts . Conversion of frequencies to normal curve deviates and 

^This research was jointly supported in part by Princeton University, the 
Office of Naval Research under contract N6onr-270-20, and the National Science 
Foundation under grant NSF G-642, and in part by the Educational Testing Service . 



assignment of weights is accomplished on a gang punch. The main procedure 
is an iterative process with repetitive cycles of computation using a 602A 
Calculating Punch and an accounting machine . 

THE GENERAL LEAST SQUARES SOLUTION TO SUCCESSIVE INTERVALS 

The experimental procedure for the method of successive intervals requires 

n stimuli to be sorted N times into (k + 1) categories on some attribute continuum. 

From this sorting procedure the number of times, f , that the ith stimulus 

ig 

was placed in the gth category can be readily obtained. The frequencies in each 

row of the n x (k + 1) table generated by these data are then cumulated so that 

each entry, F , now represents the number of times the ith stimulus appeared 
ig 

below the gth category boundary, t • These cumulated frequencies are converted 

O 

into proportions, p , and then to normal deviate values, z. .At this point, 
ig ig 

some set of weights, w. , can be applied to the normal deviates to take account 

of differences in the reliability of the various proportions . 

2 
The general least squares solution to successive intervals provides 

formulae for the n scale values, m., the n discriminal dispersions, s , and the 

1 i 

k category boundaries, t , as follows: 
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The subscript cc was introduced into the above equations to indicate the 

various cycles of approximation in an iterative procedure. Thus, if some 

tentative first estimates, t , , of the category boundaries were available to 

start the iteration, equation (1) could be solved to obtain first estimates, 

s-j , of the discriminal dispersions; equation (2) could be used to find first 

estimates, m., , of the scale values; and equations (6) and (7) could be solved 

for second estimates, t , of the category boundaries, etc. The procedure could 

g2 

be iterated until two successive t- estimates were as similar as desired, 

i.e., until [t . _, , v - t Iwas negligible . It should be noted that the coefficients 
L giOL + l) gj ^ ^ 

a., b., and c. are written solely in terms of the data and do not involve iteration 
subscripts OT ; they need to be computed only once for the entire iterative 



procedure . 

A suitable first estimate of t with which to start the iterative procedure 

g 

must now be determined. In order to be consistent with the above equations, an 
origin and unit for the t- scale should be defined so that 

k n _ 
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This definition is completely arbitrary, the successive intervals scale 

being determined only within a linear transformation. A t-scale meeting such 

requirements may be obtained as follows . If v is a set of k numbers to be 

used as possible estimates of t , then 

g 
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A set of equally spaced numbers, such as the integers from 1 to k converted 
according to equation (10), would be a convenient first estimate of the category 
boundaries . The rate of convergence may possibly be increased by doubling or 
tripling the difference between successive t-estimates, i.e., instead of using 
t in the (oC + 1) iterative cycle use (t = t + 2(t - t )• 

g(cc + 1) g (a + 1) g g(a: + 1) g 

The weights, w. , appearing in the above equations may be chosen in 

any fashion as long as w =0 and w z = when p= or p = 1 . It would 

ig ig ig 



seem reasonable, however, to select these weights so as to take account of 
differences in the reliability of various proportions . See reference 2 for a more 
detailed discussion of the choice of weights. 

TRANSFORMATIONS TO SIMPLIFY THE PUNCHED CARD PROCEDURE 

Before presenting the punched card routine, the introduction of some linear 

transformations will greatly simplify the procedure . Because of the conversion 

of equation (10), the t values in the above solution are both positive and negative, 

gOC 

as are the normal deviate values, z . Since it would be advantageous to have 

ig 
all these quantities converted to positive values, the above solution was modified 

according to the following functions: 

Z := ez + f (11) 

ig ig 

gOC gcC ' ^ ' 

where e and f are transformation constcuits chosen to make z and t all 

ig gcC 

positive and to expand the range . 

The introduction of these conversions does not affect the formula for estimates 

of the dlscriminal dispersions at all . Its form is merely changed to 

k. k 

s,^ = A.Iw. T Z. -BEw. T (13) 
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where A. and B equal a and b , respectively, withZ. substituted for z . 
1 i i i ' ig ig 

However, the transformations do produce some changes in the formula 

for estimating scale values . The m values may now be obtained in terms of 

ice 

the transformed scores as follows: 

1 
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Again, B,- and C. equal b. and c. , respectively, with Z substituted for z , 
^111 ig ig 

In summary, if the transformations given in equations (11) and (12) are 

used, equation (15) may be applied to obtain some values, P. ,which can 

be converted by equation (14) into the scale values, m . New estimates of 

i 

the t „ scale can then be obtained as follows: 

ga 
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THE PUNCHED CARD PROCEDURE 
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1 . The starting point for successive intervals analysis is an n )< k table 

of the number of times, F , that the ith stimulus was placed below the 

ig 

gth category boundary, t . Since the data are rarely available in this 

g 
form, individual scores must be tallied. If a card is made up for each of 

N individuals, with n columns representing n stimuli, the category 

(1 to k + 1) into which each stimulus was placed by an individual can 

be punched in the columns appropriate for that stimulus . It may be 

necessary to use more than one column to designate each stimulus; and 

if there are too many stimuli for one card per individual, two or more 

cards might be punched for each person. Then the cumulated frequencies, 

F. , with which each stimulus was placed below each category boundary 

may be directly tallied on the IBM 101 Electronic Statistical Machine. 

These operations may also be performed on an accotmting machine . 
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The cards used for tallying the initial scores will not be needed in the rest 

of the analysis. The cumulated frequencies, F. , are placed in detail 

cards - there being one such card for each entry in a stimulus -by-category 

table. There are actually k + 1 categories; however, since F. would 

always equal N for the last category, it need not be punched. Also, there 

will usually be missing entries in the F table, so the total number of 

ig 
detail cards will generally be less than n X k. The detail cards will be 

operated upon from this point on. In addition to F. , each card should 

also contain punches for stiraulus and category designations . 

2. A deck of conversion master cards must now be set up which will convert 

cumulated frequencies, F , into corresponding transformed normal 

ig 

deviates, Z . Accordingly, each master card should contain F. and 
ig ig 

the corresponding Z. value. It will be advantageous for later operations 

if the conversion cards also contain for each F the corresponding wei^t, 

2 '^ 

w- , and the products w Z and w Z . The conversion deck should 

ig ig ig ig ig 

also be set up so that cumulated frequencies for which the corresponding 

w is zero would be converted into X punches in the Z field. The X 

ig ig 

punch is used to indicate a conversion to zero, so that Z values with 

zero weights may be easily sorted out at a later stage of the procedure . 

This single conversion deck, then, includes the weighting system, the 

conversion to normal deviates, and the linear transformation of equation 

(11). 

3. The detail cards are now sorted into ascending order of F. , the conversion 

ig 

master cards are merged in front controlling on F. , and the quantities Z. , 

ig •■■& 

2 
w. , w Z , and w. Z are transferred from the master cards by 
ig ig ig ig ig 
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gang punching. Each detail card now contains entries of F , Z. 

ig ^S 

(where the Z. value may be an X punch), w. , w Z , w. Z , and 
ig ^ " ig ig ig ig ig 

punches identifying stimulus, i, and category, g. 

4. Some of the detail, cards will not be used in the analysis, since they contain 

little or no information . At this point, all cards with an X pimched in the 

Z. field can be removed. The X punch indicates that the information in 
ig 

the card was based upon a frequency so unreliable that it had been 

weighted zero . In case there is only one Z for any stimulus i, the card 

ig 

for this Z. is to be removed and the stimulus dropped from the study. 

5 . The following steps may now be performed on a desk calculator: 

(a) Using the integers 1 to k as v , solve equation (10) for initial 
estimates of the category boundaries, t . If the number of stimuli 

is so large that it is difficult to find the k values of ^w by 

1 ig 
n ^ 

hand, the cards may be sorted by category and Ew obtained for 

i ig 
each g on an accounting machine , 

(b) Convert the t values into positive scores, T , by equation (12). 

&1 gl 

6 . Gang punch T values into detail cards with the corresponding category 

O 

designation, g. 
7 = Sort the cards in order of stimvili and insert trailer cards . 

8. Punching a trailer card for each stimulus, run off the sums 
k k k 2 

e ig ' ff ig ig» ^"*^ ^ig ig °"^ ^" accounting machine . 

9 , Using a desk calculator, compute the coefficients A. , B. , and C- , 

according to equations (3), (4), and (5), respectively, with Z. 

substituted for z . For large numbers of stimuli these values may be 

ig 
more efficiently computed on the 602A Calculating Punch . 
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10 . Wi± the cards in order of stimuli and punching a trailer card for each 

k k 

i, compute on the 602A the values Ew T , and Hw Z T . 

g ig gl g ig ig gl 

11. Using a desk calculator, solve equations (13) and (15) for the n values of 
s. and P. , respectively. Again the 602Ais probably more efficient 
for solving these equations when the number of stimuli involved is very 
large . 

12. Punch into each detail card the corresponding s , and P., value. 

il il ^ 

13. Sort the cards in order of category, and for each g obtain Zw , 

n n i ig 

Zw. Z. s., , and Z w. P on the 602 A. 
i ig ig il ' i ig il 

14. Again at this point, hand computation is probably better for solving 
equations (16) and (17) for a new set of t values . 

15. Convert t to positive scores according to equation (12) and gang 

g2 
pvmch the T values into detail cards with corresponding category 
g2 

designations . 

16. Beginning with step 10, repeat the above computations to obtain new 

estimates s „ , P.„, and t . Steps 10 to 15 can then be iterated, 
12 12' g3 

replacing the subscript 1 with an CC appropriate to the iterative cycle 
being performed, until two successive t-estimates are as similar as 
desired. 

17. When convergence has been obtained, the final scale values, m. , can 
be computed from equation (14). 
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A SIMPLIFIED METHOD FOR THE COMPUTATION OF BISERIAL CORRELATION 
COEFFICIENTS ON THE 604 ELECTRONIC CALCULATING PUNCH^ 

K. Warner Schaie , Allan Katcher, S. Frank Miyamoto and Laura I. Crowell 

University of Washington 

A number of schemes for obtaining biserial correlation coefficients by means 
of punched card computing methods have been reported in the psychological 
and computing literatures . All of these procedures ( a representative example 
is the one proposed by Johnson*) yield only the components necessary for 
computing the final coefficient and require additional manual computation on a 
desk calculator . 

The greatest difficulty in programming a straightforward procedure lies in 
the fact that it is necessary to obtain values from a table of the normal 
probability curve in the course of the computation . While this problem can be 
readily taken care of on the IBM 650 or 704. the social scientist usually has 
access only to lower order equipment, which for his purposes is generally less 
expensive . 

In this article we will describe a one -board procedure for the IBM 604 
Electronic Calculating Punch, which will automatically perform the required 
table look-up and compute the final r. . . We selected the 604 because this 
instrument is rapidly becoming standard equipment for most basic IBM installations 
and is thus readily available . In presenting our procedure, familiarity with the 604 

and with general IBM procedures is assumed, and the terminology used agrees 

3 
with the instruction manual for the 604 . 

a 

The authors gratefully acknowledge the research support provided by a grant 
from the Agnes H. Anderson fund of the University of Washington. 

^ow with the Department of Psychiatry and Neurology, Washington University 
School of Medicine, St. Louis, Missouri. 
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Given the continuous and dichotomized scores, the mean, and the standard 
deviation for the total group on the continuous variable, the proposed method 
will summary punch, in one machine rim, one biserial correlation coefficient 
at a time . No prior information is required on the magnitute of p or q, and 
it will be shown that categorical data can be used for purposes of dichotomization 
without having to be repvuiched. This method is thus also useful in cases 
where the investigator originally intended to compute Pearsonian r's but, after 
inspection of his data, concluded that biserial correlations were more appropriate . 

The method to be described here was developed in connection with a 

research project focusing on self-concepts of communicative skills-'^ . Two 
specifically constructed questionnaires were administered to several hundred 
students who participated in three different studies . As one phase of the study, 
we were interested in examining relationships between biographical information 
and questionnaire scores . The computational example given in the final part 
of this paper is derived from this study . 
COMPUTATIONAL FORMULA 

For our purposes, the most economical formula is one suggested by 

2 
Garrett , since it requires a minimum number of constants to be carried 

during the calculation. We write the formula for the biserial correlation 

coefficient as follows: 

bis P t • ^ . 

o- t y 
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where M == the mean of the continuous variable for the larger proportion of 
P 

scores, regardless of whether this group represents or 1 scores . 

M — the mean of the total group on the continuous variable, 
<r = the standard deviation of the total group on the continuous variable. 

p = the proportion of the group having the largest number of entries . 

y = the ordinate of the normal curve corresponding to p . 
Of all the required components, a priori knowledge is needed only of Mj and 
o- J, since all other components, with the exception of y, will be obtained 
during the computation. The values of y, for a sufficient range of p's, were 
obtained from a conventional table of the normal curve . These values were then 
punched on a special deck of cards, described below. Values for M and tr 
may be computed by a number of machine routines . We used a convenient procedure 
to obtain these values on the 604, which Lunneborg, Wright, and Ax recently 
reported . 
PREPARATION OF PUNCHED CARDS 

Any set of detail cards containing scores for the continuous and dichotomous 
variables may be utilized, but a standard format will be given to correspond to 
the wiring diagrams shown in this paper . 

ColunMis 1 to 4 are assigned the subject identification code. Colunm 6 
receives an X punch, which identifies the card as a detail card and impulses 
the programs steps required for this card. The continuous scores are punched 
in columns 8, 9, and 10, In our design, three-digit scores are used for the 
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continuous variable; two-digit scores would have to be punched in the form xxO 
to conform to this layout . A zero or one is pimched for each dichotomized 
variable, one column per variable, beginning with column 37 . Where scores 
may have more than two categories, numerical scores are punched; and after 
a decision has been made as to where to dichotomize, allowance can be made 
in the punch panel wiring. 

Four tj^pes of function cards are required . The first is a card which 
precedes the deck during calculation and clears the storage units . This card 
receives an X pxmch in column 78 . 

The second function card reads in the required constants and impulses part 
of the computation. This card receives the reciprocal of the total number of 
cases (1/N) in columns 1 to 5. The mean for the total number of cases M is 
punched in colunms 8 to 11, and the standard deviation for the total group 
(T is assigned to colunms 12 to 15. Both mean and standard deviation are 
taken to two decimal places, while the reciprocal is taken to five decimals . An 
X is placed in column 79 to impulse the program steps wired to occur when 
this card passes through the calculator. 

Next we prepare the card which receives the final result. It requires a 12 

punch in column 77 . (A 12 rather than an X punch must be used to permit 

immediate transfer to avoid normal punch suppression. See diagrams.) 

This card is also prepmiched with the variable identification number to ensure 
easier recording of the results . 
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Finally we prepare the table look-up deck of 40 cards, each of which 
receives an X punch in column 76 . Each of the X76 cards further receives 
a numerical identification in columns 1 and 2 . The values for p are punched 
in columns 21 and 22, and the corresponding values for y (staggered one card 
to permit delayed pickup) are punched to three decimals in columns 23 to 25. 
Table 1 gives the required values for punching this deck, which covers splits 
from 50-50 to 90-10. This should be adequate for most empirically derived 
distributions . 
COMPUTATIONAL PROCEDURE 

In developing the present procedure, advantage was tciken of the limited 
decision- making capacity of the 604 and also of the possible use of this machine 
as a miniature card-programmed calculator . Table 2 lists all the steps required 
in the computation . These will be summarized in the following paragraphs . 

The cards are arranged in the following sequence: X78, detail cards 
(X6), X79, table look-up deck (X76), 12-77. When the X78 card is read, the 
storage units are cleared; but no calculation is otherwise permitted to occur 
on this card. 

As the first X6 (detail) card is read, a 1 is emitted into the counter . The 
dichotomized score is read through a digit selector and, depending upon 
whether the score is zero or unity, the 1 in the counter is assigned to different 
storage compartments . Depending upon this information also, the numerical 
score is read into different storage compartments . Both the coimt as well as 
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TABLE 1 

Values of p = y Punched for the Table -Look-Up Operation on the 604^ 

Column 1-2 21-22 23-25 

1 52 399 

2 53 398 

3 54 398 

4 55 397 

5 56 396 

6 57 394 

7 58 393 

8 59 391 

9 60 389 

10 61 386 

11 62 384 

12 63 381 

13 64 378 

14 65 374 

15 66 370 

16 67 366 

17 68 362 

18 69 358 

19 70 353 

20 71 348 

21 72 342 

22 73 337 

23 74 331 

24 75 324 

25 76 318 

26 77 311 

27 78 304 

28 79 296 

29 80 288 

30 81 280 

31 82 271 

32 83 262 

33 84 253 

34 85 243 

35 86 233 

36 87 223 

37 88 212 

38 89 200 

39 90 188 

40 91 176 



a 



20 Staggered to permit delayed pickup 



TABLE 2 

Program Chart for the Biserial Correlation Computing Procedure 



Card 


Prog. 








Units 




No. 


Step 


Read Out 


Read In 


Function 


Into Out of 


Suppression 


6 


Read 




Mult Quot 








79 


Read 




Mult Quot, 
FS2, FS4 








76 


Read 




FSl, GS2 








78 


1 


Counter 


FSl, FS3 


Positive 




NX78 


78 


2 


Counter 


GSl-2, GS3-4 


Positive 




NX78 


6 


3 




Counter 


Emit 1 Positive 


NX6 


6 


4 


FSl or FSS'^ 


Counter 


Positive 




NX6 


6 


5 


& Reset Counter 


FSl or FS3^ 






NX6 


6 


6 


Mult Quot 


Counter 


Positive 




NX6 


6 


7 


GSl-2orGS3-4^ 


Counter 


Positive 




NX6 


6 


8 


Counter 


GSl-2 or GS3- 


4^ 




NX6 


79 


9 


FSl 


Counter 


Positive 




NX79 


79 


10 


FS3 


Counter 


Minus, Balance 


NX79 










test for step 














suppression 






79 


11 


& Reset Coimter 








NX79 


79 


12 


FSl 




Multiply Positive 


NX79 & 
minus 


79 


13 


FS3 




Multiply Positive 


NX79 & 














plus 


79 


14 






1/2 adjust 


3 


NX79 


79 


15 


FSl 


FS3 






NX79 & 
minus 


79 


16 


& Reset Counter 


FSl 




4 


NX79 


79 


17 


GSl-2 


Counter 


Positive 


3 


NX79 & 
minus 


79 


18 


GS3-4 


Counter 


Positive 


3 


NX79& 
plus 


79 


19 


FS3 




Divide 




NX79 


79 


20 


& Reset Counter 








NX79 


79 


21 


Mult Quot 


Counter 


Positive 




NX79 


79 


22 


FS2 


Counter 


Negative 




NX79 


79 


23 






1/2 Adjust 




NX79 


79 


24 


& Reset Counter 


GSl-2 


Positive 


2 


NX79 


79 


25 


GSl-2 


Counter 


Positive 


5 


NX79 


79 


26 


FS4 




Divide 




NX79 


79 


27 


& Reset Counter 








NX79 


79 


28 


Mult Quot 


GS3-4 






NX79 


79 


29 


FSl 


FS3 






NX79 


79 


30 


GS3-4 


Counter 


Negative 




NX79 & 
minus 


79 


31 


& Reset Counter 


GS3-4 






NX79& 
minus 


76 


32 


FS3 


Counter 


Positive 




NX76 


76 


33 


FSl 


Counter 


Negative 




NX76 


76 


34 






Zero Check 




NX76 


76 


35 


GS2 


FS4 






NX761' 


76 


36 


& Reset Counter 








NX76 


77 


37 


FS3 


Counter 


Positive 


5 


NX77 


77 


38 


FS4 




Divide 




NX77 


77 


39 


& Reset Counter 








NX77 


77 


40 


GS4 




Multiply Positive 


NX77 


77 


41 






1/2 Adjust 


2 


NX77 


77 


42 


& Reset Counter 


GS3-4 




3 


NX77 


77 


43 


GS3-4 


Coimter 


Positive 




NX77 


77 


Punch 


& Reset Counter 











^The dichotomized score read in through a digit selector determines the position of calculator 
selector 1. If a zero or low-order digit is read, FSl and GSl-2 are used; if a unity or high-order 
digit is read, FS3 and GS3-4 are used (see text). 

Step 35 is suppressed also whenever a zero check impulse transfers calculator selector 2, 
through which this program is wired. 
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the numerical score is then added to the totals accumulating in the correct 

storage units. This selection is performed by impulsing calculator selector 1 

to transfer whenever a zero is read. If categorical scores have been punched in 

the detail card, it is necessary only to wire the digits which are to be considered 

zero through bus hubs to the calculator selector pickup to achieve the same 

result as if a zero had actually been punched in the card. 

The X79 reads in the information on M , o- , and 1/N,as well as impulsing 

a number of logical and arithmetical steps . The total counts for zero and unity 

entries are compared. The larger count is then divided by its N to obtain 

M , from which M is next subtracted. The remainder is then divided by 
P t 

a and stored. Depending on the outcome of the balance test, an adjustment 
is made to give this remainder the correct sign. 

On each X76 card the p entry stored in the calculator is compared with 
the one read from the individual X76 card. A zero check impulse is emitted 
whenever the two do not compare, thus transferring calculator selector 2 
to inhibit reading in the y value. At the time when the two amounts do agree, 
the corresponding y value is permitted to read into calculator storage from 
the X76 card. Since the zero check impulse is emitted one cycle later, the y 
entries are staggered by one card on the table look-up deck . 

The X77 card permits completion of the calculation, p is divided by y 
and the result is multiplied by the stored product of the first part of the 
equation which was computed while the X79 card passed through the calculator. 
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The result punches on this card in columns 71 to 74, the coefficient being 
punched to four decimals . 

Figure 1 shows the punch panel wiring for this operation. It should be noted 
that all wiring is permement except for the wire from the common hub of pilot 
selector 2. This wire is moved after each rvm if more than one dichotomous 
score is punched on the detail cards, as will ordinarily be the case . Figure 2 
gives the calculator panel wiring, except for the necessary program suppressions , 
These are shown separately in figure 3 for greater clarity. 

The time required to compute one coefficient is determined by the number 
of detail cards (N) plus the required 43 function cards . It can be determined by 
the equation 

Machine time for one r = ^ '*" 3^ minutes 

bis IQQ 

Because three-digit storage units are used to record the unit count, the upper 
limit for N is equal to p = 999 or, in general, to approximately 1200 cases, 
depending upon the extremity of the splits . The method will prove most useful 
where there are at least several hundred subjects and where there are few 
continuous and many dichotomous variables . 
COMPUTATIONAL EXAMPLE 

As already mentioned, the above method was developed for a study of self- 
concepts of communication skills . As a simple illustration we shall take the 
scores for six subjects on the commvmication scale and the dichotomized score 
on "number of siblings . " The communication scale required the subject to 
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Figure 2. Calculator Panel for Biserial Correlation 
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Figure 3. Program Suppression for Biserial Correlation 



appraise his performance as a communicator in a wide variety of situations . 

A scoring scheme was employed so that total scores on the scale reflected 

the extent to which a subject felt his over -all performance was high or low. 

Table 3 gives the information for these subjects punched into the detail cards 

and the X79 card. 

TABLE 3 
Detail Cards for the Computational Example 



Column 



Subject 


Communication 


No. of 


Identif. 


Score 


Siblings 


1-4 


8-10 


37 


0001 


075 




0002 


085 




0003 


110 




0004 


080 




0005 


090 


2 


0006 


125 


3 



X =94.17 
o- =17.65 
1/N =.16667 

The categorical variable (number of siblings) was dichotomized into subjects 
with one sibling and subjects with more than one sibling. Since a punch of 1 in 
the dichotomized variable column was the low -order digit, the 1 on the digit 
emitter was wired instead of zero . In all other respects the panels remained 
as described. The following paragraphs describe in detail the program steps 
of the required computations . 

Only program steps 3 through 8 are permitted to operate on the detail 
cards. On program 3, a 1 is emitted into the counter. Since a low-order digit 
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is read on card 0001, the 1 in the counter is transferred to FSl on program 5. 
(Nothing is added to the counter on program 4, since FSl has been cleared by 
the X78 card preceding the first detail card.) Similarly, the continuous score 
(75) is transferred on program 8 to GSl-2. On the second detail card (since 
again a low-order digit is read), program 4 adds the stored count of 1 to the 
counter and program 5 returns the augmented count of 2 to FSl . Similarly, 
program 7 adds the accumulated sum of the continuous scores to the one read 
in on card 0002 and program 8 returns the augmented sum to storage . The same 
process is repeated for cards 0003 and 0004, as both have low-order digits on 
the dichotomized variable. A high-order digit is read on card 0005, and therefore 
no impulse is available to transfer calculator selector 1 . As a result, the unit 
count is stored in FS3 instead of FSl, and the sum of continuous scores is 
accumulated in GS3-4 instead of GSl-2. The information from card 0006 is treated 
in the same manner. 

Table 4 gives the values which appear in the storage units and counter 
after completion of each program step, beginning with step 8 . These program 
steps will be discussed in detail below. 

The X79 card has the values given in table 3. Of these, M is stored in 
FS2, <r in FS4, and 1/N in MQ. 

On program 9 the sum of the low-order digit count enters into the counter 
from FSl . On program 10 the sum of the high-order digit count enters 
negatively from FS3 and a program test for balance suppression is performed. 
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TABLE 4 

Computational Example 



Step 


Factor Storage 


MQ 


Counter 


General Storage 


1 


2 


3 


4 


1 2 


3 4 


Computations During the X79 Card 


8 


4 




2 






4 


350 


215 


9 




94.17 




17.65 


.16667 


4 






10 












2 






11 



















12 












.66668 






13 


















14 












.67168 






15 






4 












16 


67 

















17 












350000 






18 


















19 










87.500 









20 


















21 












87.500 






22 












-6.670 






23 












-6.675 






24 















-6.67 




25 












-6.67000 






26 










-.376 


Rem. 






27 



















28 
















-.376 


29 






67 












30 












.376 






31 

















.376 


Computations During the 17th X76 Card 


32 


68 










67 


.362 




33 












-.01 






34 












-.02 






35 










.362 








36 



















Computations During the 77 Card 


37 


91 










.670000 


.176 




38 










1.851 


•Rem. 






39 



















40 












.696776 






41 












.696826 






42 

















.6968 


43 












.6968 







All values are shown when they first enter the counter or storage units, remaining 
there imtil a change is indicated. 
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The result of this test determines which of several alternate program steps 
will be suppressed (see table 2). In our example, a plus balance results 
and program levels 13 and 18 are suppressed. Program 11 clears the counter; 
and on 12 the count for the larger sub-group is multiplied by the reciprocal of 
the total number of cases to yield p, the proportion of the larger group. This 
result is rounded off on program 14 and, on program 16, stored in FS2 to 
two decimal places. Program 15 serves to shift the N for the larger sub-group 
from FSl to FS3 and is required to make available additional storage facilities. 
On program 17 the sum of the continuous scores for p is read into the counter 
and is divided, on program 19.. by the N for this group to yield the mean for 
the continuous scores for the larger sub-group . After the result has been 
transferred back into the counter, M is subtracted by negative entry from 
FS2 , The result of this operation is rounded off in the last place and then stored 
to two decimal places in GSl-2. On program 26 it is read out again into the 
counter and divided by o- . The quotient is stored in GS3-4; but its sign is 
reversed on programs 30 and 31, since the low -order digit group is the larger, 
to assure that the final coefficient will have the correct sign. The program levels 
passed over in the description serve to clear and shift the required storage 
units . 

The next step in the calculation consists of reading in the corresponding 
value of y for the value of p which has been computed while the X79 card 
passed through the calculator. The deck of X76 cards constitutes the table 
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required for this purpose . Successive trial values of p are read from the X76 
cards and are compared with the computed value of p stored in the calculator. 
Whenever the two values balance, the corresponding value for y is stored in 
the calculator . On step 32 the computed value of p is read into the counter from 
FS3. Then on 33 the read-in trial value is subtracted. On 34 a zero check 
test IS performed. 

We found the value of p for our example to be . 67 . Therefore, no changes 
occur in the calculator mitil card # 16 of the X76 deck passes. No zero check 
impulse is emitted because observed and trial values agree; and on the next 
card the y value is entered into FS4, since step 35 is permitted to be active . 
Table 4 shows the values found after calculation is complete on X76 #17. 

Since a zero check impulse will be emitted on all subsequent X76 cards, 
no further changes occur until the 77 card enters the calculator and impulses 
the remaining computations . On step 37 p is read into the counter into fifth 
place and is divided by y on step 38 . After resetting the counter this result 
is multiplied on step 40 by the product of the preceding calculations stored in 
GS3-4. This final product is then rounded in the second position on step 41. 
The last two steps reduce the coefficient to four decimals, which are punched on 
the 77 card to give the biserial correlation coefficient of . 6968 for our example . 
SUMMARY 

This paper describes a method for computing biserial correlation coefficients 
on the 604 Electronic Calculating Punch in a single -step procedure, provided 
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the mean and standard deviation for the total sample are known for the continuous 
score . The method is also adaptable to the case where categorical scores must 
first be dichotomized before biserial r can be computed. A computational 
example as well as all required wiring diagrams and program charts are given . 
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A METHOD FOR THE PACKAGED PROCESSING OF A STATISTICAL 

ANALYSIS ON THE IBM 650 

Leon H. Somerall and Nicholas A. Habibe 

Air Weather Service 
Asheville, North Carolina 

The primary purpose of this paper is the presentation of a logical 
program design for the "one package" processing or continuous solution of 
a standard statistical problem on the IBM 650 Magnetic Drum Data 
Processing Machine . 

The one package concept embodies all the functions (except the pre- 
liminary sorting or assembly and the printing of the output) formerly attained 
through "step" or discontinuous procedures requiring the use of diverse IBM 
equipment such as sorters, summary punches, collators, accounting machines, 
and/or calculating equipment such as calculating punches, card-programmed 
calculators , Proper program logic enables the 650 to simulate functions 
featured in component- type punched-card equipment. 

The standard statistical problem to be discussed consists basically of 
two parts: (1) performance of a bivariate frequency distribution by a table 
look-up method and (2) evaluation of some statistics pertaining to the distribution. 

The printed result, figure 1, Shows two families of imposed class 
intervals X, Y, within which the single -valued variables or observations x, y 
will fall uniquely and respectively, and the family of groups Z, called the 

parameter of the distribution. Each member of the family Z contains all the 

a 
members of X and Y. 

Part 1 seeks to build, by groups of data, a drum table of frequencies 

which will reflect the fitting of each x, each y, and the pair (x, y) into the 

appropriate class interval member X., Y., and (X., Y), respectively, to add 

a count each in the three drum cells containing the corresponding frequencies , 

3. 

The terminology and notation at the end of the paper will be found helpful . ^3 
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A count is also added into a cell reserved to keep track of how many valid 

observations have been processed. Thus, the table will furnish a running 

record of how many x's, how many y's, and how many (x, y)'s fall respectively 

in each class X. for all values of Y, in each class Y for all values of X, and 

1 J 

in each pair (X. Y.). The members of X and Y may be in random sequence 

and not necessarily equi-spaced. If the class intervals are equi-spaced and 

continuous, the distribution will not necessarily require table look-up. 

Each group Z. in the printed output, figure 1, represents a corresponding 
group Zj, of input cards, shown in figure 2 arranged in groups Z , Z , . . . , Z . 
Prior to processing in the 650, the cards are assembled in groups of the 
parameter Z by means of a sorting operation. Each card contains three fields: 
group Zj^ to which the card belongs, the observation x, and the observation y. 

Part 2 is concerned with the determination of means, standard deviations, 
and correlation coefficients. Of course, any other statistics may be included 
by appropriate programming. 

The printed output, a modified replica of the drum table, is obtained by 
processing in the IBM 407 Accounting Machine the punched output of the IBM 650. 
The headings X, , X , . . . , X^ in figure 1 are not present in the drum table, 
figure 4. They are printed by 407 control panel wiring. The drum area 
assigned to the table of frequencies and statistics must contain (m+3) (n+1) 
cells to house the frequencies, the identification, and the statistics . On the 
printed output, X runs horizontally and Y runs vertically. On the drum table, 
however, the X and Y directions are interchanged to facilitate 650 punching 
and 407 printing of the output one line per card . 

The program design is based on the general scheme of effecting the 
following data analysis on each valid card of any group Zj^: 
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(a) Validity check on x, y 

(b) Relating each member X. of X to the family Y by means of the 

frequency f., , which is the sum of the frequencies f.„ , f.„, ..... f. , 
Ik ^ n ilk i2k ink 



n 



in any one class X. for all values of Y- f = Y f 

' ^ j = l ijk" 

(c) Relating each member Y. of Y to the family X by means of the frequency 



git, which is the sum of the frequencies f, ., , f„., , .... f 
J^ Ijk 2jk mjk 



m 
iilK' ' mik 

m 

any one class Y. for all values of X: g. = ]^ f . 



jk i ^ 1 ijk 

(d) Relating each member X. of X to each member Y of Y by means of 

1 J 

the frequency f^^j, of the simultaneous pair (X., Y.). 

Each pair (x, y) is considered valid if x and y comply with some 

imposed conditions of restraint. An invalid card is disregarded. 

The frequencies or counts f., , g., , and f. ., are properly recorded 

on the drum table as a result of an interesting table look-up technique 

to be explained below. Thus, table look-up is the principal data 

processing link between the information contained in the input cards 

and the corresponding counts to be recorded in the drum table of 

frequencies . 

2 2 

(e) Progressive recording of X^x, Yy> Z^ , Zy 'Z^y, N for the evaluation, 

by groups, of the mean (5L ) of all x's, the mean (y, ) of all y's, the 

standard deviation (o- , ) of all x's, the standard deviation (o- ) of 

xk yk 

all y's, and the correlation coefficient r. between all the x's and all 
the y's in each group. 
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Figure 2. Batch or Deck-Sorted or Assembled in Groups Z^ Z , . . ., Z 

PROGRAM LOGIC 

An overall procedure for the logical processing of data is shown in the 
block diagram of figure 3, which conforms to the block diagramming conventions 
outlined in the manual "IBM 650 Problem -Planning Aids . " It is generally 
applicable in situations where it is desired to cause the production of a set of 
outputs whenever a difference exists between the Z value of two adjacent cards . 
This output will provide the result of the analysis of a Z group . 

The following discussion ties in the block diagram, figure 3, and the 
drum display, figure 4. 
BLOCKS 

01 . Early during the conception of the flow chart, a distinct fact is 
recognized: separate reading instructions must be provided for the first card 
only; the reading of any other card but the first must be accomplished with a 
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different set of instructions, as provided by block 05. Exclusive reading 

instructions for the first cards is dictated by (a) the uniqueness of the position 

of this card in the deck, (b) the necessity of clearing (block 02) the drum area 

assigned to the table of frequencies before the data analysis is to begin, and 

(c) the prevention of an output routine at the start of the data process . An 

instruction in block 01 reads the first card of the first group onto the drum, 

i.e., the observations x, y and the identification Z enter the read input 

area where the Z is then called Z . 

r 

02 . Storage clearing is a prerequisite before the analysis of each group 
Z, starts . Instructions in this block clear to zero all the memory locations 
assigned to cumulative data (frequencies or counts and other summations) . It 
is generally faster to clear storage during the output routine . If the amount 
of time devoted to output computations of statistics in block 07 is less than 

the time available for maximum punch output rate, then storage clearing can be 
incorporated in the output routine to take advantage of the remaining available 
time. In this case, block 02 is operative for the first card of the first group 
only instead of the first card of each group. In other words, when storage 
clearing takes place at output time, block 02 is still operative for the very 
first card because there is no convenient way of divorcing storage clear in 
the output routine to have it apply to the very first card. 

03. This set of instructions stores the first Z , i.e., the Z of the first 

r 

card read in each group, in another location L(Z ) where it is renamed Z . 

S o 

The purpose of the transfer is two-fold: (1) to compare the Z of the first 
card of every group with the Z's of subsequent cards in order to detect a 
change in Z or transition to a new group, and (2) to identify the output. The 
use of the punch output area for the storage of Z permits most conveniently 
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the performance of the two -fold function of Z . 

s 

04. In this particular problem, data analysis represents the following 
collection of operations applied to each card: (a) validity check, (b) progressive 
recording of cumulative frequencies by a table look-up method, (c) progressive 
recording of summations for the evaluation of statistics . Other operations, 
based on the individual values of the read elements, could be performed during 
this data analysis phase . 

A validity check, a test to determine if the read variables x, y fulfill 
certain criteria or conditions of restraint, precedes data analysis . When the 
imposed conditions are not met, the card is disregarded and the program 
logic proceeds to read another card. 

Progressive recording of cumulative frequencies for any valid card is 
accomplished through the following routine: (a) at the beginning of each group 
the drum area assigned to cumulative data is cleared to zero, as explained 
above; (b) a table look-up operation selects the appropriate terminal value 
(X(.)j of the class interval X^ in which the read x belongs, to yield the address 
L(fijj) of fjjj, that is, one of the m locations assigned for the accumulation of 
the frequencies in the class X. for all values of Y; the address being known, 
f ., is sent to the accumulator to be increased by a coimt of one and then returned 
to its original location; (c) a similar routine for the read y yields the address 
of g , that is, one of the n locations assigned for the accumulation of the 
frequencies in the class Y. for all values of X; the contents of this address are 
then increased by a count of one as explained above; (d) the address of f. ., , 

IJK 

f, that is, one of the mn locations for the accumulation of the frequencies in 



r 



the pair (X., Y.), is attained via the addresses of f., and g , ; f is then 
1 J Ik ^jk ijk 

increased by a count of one . 
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2 2 

Progressive recording of ^x, ^y, ^x , J]y , Y,^> ^^"^ ^ "^^^^ ^°^ 

involve table look-up because these summations have a fixed reference (or 
data address) throughout the data process . The search for addresses mentioned 
in block 04 (b), (c), and (d) requires a variable address scheme capable of 
giving all the locations of the table of frequencies as a function of the read 
variables . The instructions to store the counts or frequencies are automatically 
modified for each card in accordance with the individual values of x, y. A 
number of schemes could be devised to accomplish address searching, some 
of which may or may not involve table look-up and/or computation. In the pre- 
sent problem, TLU and computation are featured in tne determination of 
variable addresses because of the flexibility afforded by such a scheme. 

In the conventional sense, TLU is the process of best matching a given 
argument to one of a collection of tabulated arguments and reading off the 
tabulated function associated with the tabulated argument . In the 650 sense, 
however, the standard TLU operation yields the address or location of one 
out of many tab arguments equal to or next larger (if no equal exists) than the 
given argument. The tab function, stored with or a constant number of locations 
away from its corresponding tab argument, is then extracted by making use 
of the found tab argument address . The tab function could be any function, 
such as a function of the argument, an instruction, an address, a constant, etc. 

The standard TLU operation is susceptible to many refinements to 
accomplish address searching more efficiently. In the TLU teclmique presented 
in this paper, the tab function can be computed as a function of the location of 
the corresponding tab argument. Then, because no extraction takes place, 
drum area for the storage of tab functions is unnecessary. Since in this problem 
there are two given arguments x and y, two families of tab arguments will be 
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required. The families of terminal values X and Y of the class intervals X 

t t 

and Y, which are the required tab arguments, are referred to here as the 
X directory and the Y directory. The location of the frequencies in the drum 
table of frequencies constitute the "tab functions . " 

The following drum areas are required for the execution of the table look- 
up techniques: 

Contents 



Area Name 



Read input 



X directory 



Y directory 



TLU constants 



A value of x 
A value of y 
Terminal values 

Terminal values 
(Y) , (Y ),..., (Y) 

t 1 t 2 tn 

^1' ^2' ^3 



Role 

Given x argument 
Given y argument 



Tab X arguments 



Tab Y arguments 

Computation of 
W^). Ug.^), L(f..^) 



In order to store counts for the frequency in each class X. for all values 

of Y, in each class Y. for all values of X, and in each pair (X., Y.), the 

locations L(f., ), L(g. ), and L(i. ) must be computed for every card with 
iic jk ijK 

valid X and y using the equations 

L(fij,)= L(Xp. -f k^ 



L(gjj^)=50L(Y^). + k^ 



(1) 
(2) 

L<fijk>= ^(fik> + L(g.p + k3 (3) 

The numerical values of k, , k , k- are dependent on the relative position 

•"^ 2 "J 

of the drum areas assigned to the directories and to the table of frequencies . 

The sequence of events for obtaining the unknown addresses in equations (1), 
(2), (3) and for storing the counts in each address follows: 
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Step 


Operation 


Result 


1. 


TLU 


Read input x vs . X directory 


MX), 


2. 


TLU 


Read input y vs . Y directory 


L(Y). 


3. 


Computation 


Equation (1) 


1-<W 


4. 


Computation 


Equation (2) 


Ugjk> 


5. 


Computation 


Equation (3) 


l-tfijk) 


6. 


Store count 


One unit added to contents of 


'ik 


7. 


Store count 


One unit added to contents of 


^jk 






i-'SjiP 


8. 


Store count 


One unit added to contents of 
L(f . ) 


^ijk 



To illustrate the table look-up technique presented, suppose 

(a) a card in the first group (k = l) undergoing data analysis contains the 
following: x = x , y = y , Z = Z , x and y valid. 

(b) X falls in the class interval X (i = 5), whose terminal value is (X ) ; 

-* 5 t 5 

y falls in the class interval Y , (j = 6), whose terminal value is (Y ) - . 
o o to 

(c) the family X contains nine class intervals, the family Y contains 
seven, and the drum area layout is the one shown in figure 5 . 

The problem is 

(a) find the drum location L(f ) of the frequency f , i.e., the location 

Di 51 

(0768) of the frequencies of occurrence in the class X , group Z , 

5 1 

regardless of Y. 

(b) find the location L(g^ ) of the frequency g. . i. e., the location (0673) 

61 61 

of the frequencies of occurrence in the class Y^, group Z , regardless 
of X. 

(c) find the location of L(f ) of the frequency f_^ , i.e., the location 

561 561 
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Constants 



X 

Y 



"3^ 



-22,077 



-773 



X. directory 



/ 



0400 

(xp, 


0450 

(^t>i 


0401 
(Xp2 


0451 


0402 
(XP3 


0452 


0403 
(XP4 


0453 


0404 
(Xt)5 


0454 


0405 

(xpe 


0455 


0406 


0456 
(Yj), 


0407 




0408 





^Y, directory 



Y- 


— >- 
















Z,,Y, 


2r^2 


ZpYg 


Zl,Y^ 


ZpYs 


^I'Ye 


ZpY^ 


h 


h 


0414 
^111 


0464 
'l21 


0514 
^131 


0564 
'l41 


0614 
^151 


0664 
^161 


0714 
^171 


0764 


^1 


^211 


f 

^221 


^231 


• • • 


• ■ • 






'21 


h 


^311 




• • * 










'31 


h 
















'41 


"i 












0668 
^561 




0768 
'51 


h 
















'61 


















'71 


















'81 


















'91 




0423 


0473 
^21 


0523 
^31 


0573 
g41 


0623 
^51 


0673 
^61 


0723 

8,1 


N 





Figure 5. Drum Layout for the Illustrative Example 
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(0668) of the frequency of occurrence in the pair of classes (X Y ), 

5 6 

(d) increase by one count each, the contents of the locations L(f ), 

L(g, ), and L(f ) . 
^61" ^ 561 

The solution is: 

Step 1 , TLU of given argument x^ versus the X directory yields 0404, the 

o t 

location of the tab argument (X ) equal to or next higher than x . 

Step 2 . TLU of given argument y , versus the Y directory yields 0455, 

6 t 

the location of the tab argument (Y. ) equal to or next higher than 

t 6 

ye- 

Step 3. Equation (1) gives: 

L(f ) = 0404 4- 364 = 0768 
51 

Step 4 . Equation (2) gives: 

L(g^ , ) = 50(0455) - 22, 077 = 0673 
61 

Step 5 . Equation (3) gives: 

W^(^j)= 0768 + 0673 - 773 = 0668 
Steps 6 through 8 . The contents of locations 0768, 0673, and 0668 are 

increased by one unit . 

05. The program logic shown by the block diagram requires two blocks 
of instructions for the 650 to execute card reading. Block 05 reads all cards 
except the first card of the first group, which is read by block 01. 

06 . This test compares the identification Z of the first card of a group 

s 

with the identification of each card within the group . If the equal condition 
exists, the card read (just before the test) is part of the group currently under 
process, in which case the data analysis and card reading will not be disturbed 
until an unequal condition is detected. If the test yields an unequal result, the 
card read (just before the test) is the first of a new group, Implying that the 
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data collected and analyzed on an individual card basis is complete for the group; 
card reading is then temporarily interrupted until the output routine takes place . 

07. The output routine incorporates two subroutines: (a) computation or 
further data analysis on a group basis and (b) assembly of the data in the punch 
output area. Means, standard deviations, correlation coefficients, etc, are 
computed in subroutine (a) . The data to be punched is transferred from the 
table of frequencies to a punch output area by subroutine (b) . It will generally 
require the computation of the data address (especially if a "loop" is used) 
of the instruction to "go and get" the appropriate data items for transfer to 
the output area. The data address of this instruction can very conveniently be 
inserted in a "clear to zero" instruction so that a "write and erase" subroutine 
can be developed. Upon completion of the output routine the data process starts 
all over again, this time with block 02 or 03 . 

The generalized program described above can be easily modified to do 
a variety of accounting and statistical problems . 
TERMINOLOGY AND NOTATION 

Conditions of Restraint Criteria for acceptance of data or restrictions 

imposed in the course of a process . 

Batch or Deck A collection of groups of cards . 

Group A collection of items or cards having a common 

characteristic (such as identification) . 

Valid Data which meet the criteria for acceptance 

or fulfill the imposed restrictions . 

Table Look-up The process of best matching a given argument 

with one of a collection of tabulated arguments 
to read off the tabulated function associated 
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Given Arg 



Tab Arg 

Tab Function 
Directory 



Extraction 



Class Interval 



X, y 



X 



with the tabulated argument . 

Given argument, a value of the variable which, 

by comparison with a tabulated argument, is 

used to search for the tabulated function. 

Tabulated argument, a value which best matches 

the given argument. 

The value associated with the tab argument. 

A relatively small table or drum area consisting 

of tab arg, tab functions, constants, etc., 

sometimes necessary for access to a larger 

drum area or table . 

The process of obtaining the tab function by 

means of 650 coded steps , 

A range of values comprising and contained 

within an initial value and a terminal value . 

Observations or variables to be distributed 

in class intervals. They have a finite thou^ 

variable number of digits and belong in one 

class interval only, provided the class intervals 

do not overlap. 

Family of class intervals in random sequence 

and not necessarily equi-spaced, representing 

all the members X, , X , . . . , X . 
1 2 m 
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X. 



X. 



<Vi 



a,\ 



Y. 
J 



<M 



Any one member of the family X, one class 

interval whose range is from (X„) to (X ). 

U 1 1 1 

(i=l, 2, ..., m). 

Family of initial values of the family X, 

representing all the members (X )., (X ) , (X ) . 

10 2 0m 

Any one member of the family X^., the Initial 
value of the class interval X. . 
Family of terminal values of the family X, 
representing all the members (X.) , (X )2, 

..., (X) . 

t m 

Any one member of the family X , the terminal 

value of the class interval X . 

i 

Family of class intervals in random sequence 

and not necessarily equi-spaced, representing 

all the members Y , Y , . . . , Y . 
12 n 

Any one member of the family Y, one class 

interval whose range is from (Y ). to (Y ) 

OJ tj 

(j = l, 2, , . ., n). 

Family of initial values of the family Y, 

representing all the members (Y ),, (Y ) , 

••-<v„- 

Any one member of the family Y , the initial 

value of the class interval Y. . 

J 

Family of terminal values of the family Y, 
representing all the members (Y ),, (YJ2, • • •» (Y^jj 
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(Y ). Any one member of the family Y , the terminal 

value of the class interval Y . 

J 

Z Family of group identifications representing 

all the members Z,, Z„, .. ., Z . 
12 p 

Zj^ Any one member of the family Z (k= 1, 2, . . ., p). 

f . .1^ A frequency or number of times the pair of 

observations (x, y) falls within the pair of 

class intervals (X., Y.) in a group Z, . 
n ^ 

ik = A f , any one member 

j = l ijk 

n n n 

f = y f , f = T f , ..-. f = y f 

Ik ^ lik 2k ^ 2jk rnk ^ mjk 

j-1 ' j=l j=l 

of a family f of summations of frequencies 

f . . falling within any one class interval X. for 

all values of Y in a group Z . 

m 
g^jj ~ y f ^^y °^^ member 

i=l ijk' 

mm m 

of a family g of summations of frequencies f . 

falling within any one class interval Y. for all 

values of X in a group Z . 
X, Arithmetic mean of all x's in group k. 

yj. Arithmetic mean of all y's in a group k. 

"" xk Standard deviation of all x's in a group k. 
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<T 



yk Standard deviation of all y's in group k. 



r. Correlation coefficient of y vs . x in a group k . 

Z Group identification of any card just read. 

Z„ Group identification of any group's first card 

just stored. 
L(U) Drum address or location of U, where U is 

any quantity . 
i = 1, 2, . . ., m. 

j = 1, 2, . . ., n. 

o Initial . 

r Read (past tense) . 

s Stored. 

t Terminal . 
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AN INTERPRETIVE SUBROUTINE FOR THE SOLUTION OF 

SYSTEMS OF FIRST ORDER ORDINARY DIFFERENTIAL 

EQUATIONS ON THE 650 

Franz Edelman 
RCA Laboratories 



INTRODUCTION 

The problem of solving ordinary differential equations occurs very 

frequently in technical computations . The writing and subsequent checking of 
such programs is extremely wasteful of both machine and programming time . 
Because of this, a general purpose program has been written which will 
automatically perform this function . 

To use this subroutine, the programmer need specify only initial conditions, 
the equations to be solved and their number, as well as the precision required 
in the solution. 

• The programmer has a choice of two distinct methods for solution, those 
of Runge-Kutta-Gill and Milne. The program for the latter provides automatic 
choice of optimum intervals . A choice for a fixed or floating point precision 
criterion is also provided. 

The maximum number of equations which may be solved simultaneously 
is limited only by the memory capacity of the machine and is in the neig^orhood 

The program ia written for the Interpretive System described by V. M. 
Wolontis of the Bell Telephone Laboratories, Murray Hill, New Jersey, in 
IBM Applied Science Division Technical Newsletter No.ll . 
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of thirty . 

The basic logic of this subroutine is in no way restricted to be applied 
solely to a 650 program. Therefore it is felt that, because of the flexibility and 
compactness of the logic, this program will lend itself well to transcription 
into other machine or pseudo -languages. Complete flow diagrams are furnished 
for this purpose at the end of the report. 

In designing this general purpose program an effort has been made to 
anticipate as many of the likely requirements and needs as possible, particularly 
with regard to generality, flexibility, and convenience of use . 

There are several phases of a typical differential equation problem where 
flexibility is nothing less than essential; 

The Method of Solution . In general, integration methods which keep track 
of truncation error, such as that of Milne ^, are preferable to those which do 
not, for obvious reasons . There are, however, two disadvantages to the former: 
first, a starting procedure is required before the method can be applied; and, 
second, under some conditions this method may become unstable . The first 
shortcoming is only an inconvenience, but the second renders the method useless 

For these reasons a second method is available, which is that of Runge-Kutta- 

2 
Gill . This method also doubles as the starting procedure for the Milne method. 

Output . At each step of the integration the type and amount of the output as 

well as the format are completely within the control of the user of the program . 

All the information required by the user as output is available from a given set 

of locations arranged in such a manner as not to tax his memory unnecessarily. 
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The user may, if he wishes, change the nature and format of his output during 
a run depending on some predetermined criterion. 

The Criterion for Terminating a Solution. Here again is something which 
varies from one problem to the next. Therefore the program is designed in 
such a way that the user may program this criterion, which may depend on any 
or all of the variables, to suit his purpose. 

Precision. This consideration applies only to the Milne method. The manner 
in which the precision is computed is important, since its proper choice may 
in some cases reduce machine time by more than 50 percent. This program 
offers two choices . The increasing or decreasing of the integration interval 
is governed by a comparison of the predicted and corrected values of the 
dependent variables using a predetermined number of either (1) decimal 
places or (2) significant figures . In either case this number is specified by the 
user as a floating point integer . If the solutions oscillate about zero with 
an amplitude of order unity, for example, then criterion (1) is considerably 
faster. Thus if at some point the solution is very near zero, criterion (1) would 
remain undisturbed while criterion (2), noting the poor agreement between 
predicted and corrected values (in terms of significant figures), would 
proceed to reduce the integration interval unnecessarily. On the other hand, if 
the solutions run uniformly to very large or very small numbers, criterion (2) 
is faster and more realistic. In other words, criterion (2) is always safe, even 
though sometimes excessively time-consuming. The solutions to y" + y = 
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and y" - y =0, for example, illustrate this point. 

Reducing the Integration Interval . In general, the variations of the dependent 
variables and their derivatives are not violent enough to warrant such a drastic 
reduction of the integration interval as an order of magnitude for example . 
Nevertheless it may happen that the customary reduction by the factor 0.5 
is insufficient , Therefore this factor has been left to the discretion of the user . 
Unless otherwise instructed, however, the program will use the factor 0,5. 

Increasing the Integration Interval . When the precision of the computations 
becomes too large, implying that the integration interval has become too small, 
the subroutine will automatically double the interval imtil the precision reaches 
its prescribed range. 

The above-mentioned features are believed to make this program flexible 
enough to cover a wide range of problems . 

For a better understanding of what is to follow, a list of definitions 
is given below: 

y is the independent variable; y is its (j + 1)™ value. 

y., i =1, 2, . . .,N, are the dependent variables; y.. are their values at 

f.. — (dy./dy )., 1 ■= 0, 1, 2, ...,N, are the derivatives of y. with respect 

to y„ at y„. . Thus f„ *" 1 for all j . 
■^0 ■'Oj Oj 

k.., r.., q.. are defined analytically in the following section, 
h is a positive or negative integration increment. 
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N is the number of equations to be solved. 

T is the required precision. 

H is the fraction by which h is to be reduced. 

D is the tolerance factor for doubling the integration interval . 

Ei4 is the deviation between the predicted and corrected values of y. at 
ij -^ 'i 

V„ . . E is the maximum E at y . 
OJ J ij Oj 

A table summarizing the functions of all special locations, which may be 
referred to by the user of this program, is shown at the end of this report. 

MATHEMATICAL BACKGROUND 

Any system of ordinary differential equations may be reduced, by suitable 

changes of variables, to an equivalent system of first order ordinary differential 
equations: 

y\ = \(.Yq> y^. 72' •••• y^)' <^= ^' 2. .-..n), (i) 

subject to the initial conditions 

^1^0^= ^iO' (i= i» 2,...,N). (2) 

The program described here will automatically solve such a system, 

requiring of the user only the specification of the functions f and the initial 

i 

conditions . 

For the sake of completeness the equations of both methods employed are 
given below. 
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The equations of Runge-Kutta-Gill are 



k.. = hf.. 



r. . , , — a.k. . - b.q 
1, J + 1 J iJ J ij 



^i. j+ l=yij-^ -i, j+ 1 



q- • j_ , = q- ■*■ 3r. . , - c.k.. 



>- 



j= 0, 1, 2, 3 



i=0, 1, 2, ..., N 



(3) 



(4) 



(5) 



(6) 



The initial values q.„ are taken to be zero at the first point and equal to 



q. , of the previous point, subsequently, a., b 

J J 

each given by 



a.= (i. i->4", 1+yr, 1) 



J 



= (1. i-^/ii+y^, k 



2' 3' 



c = {{, I-n/T 1+v^ i) 

The equations of Milne are the "predictor": 



(P) 4h , 



■ f + 2f 

ij i. J - 1 1' j - 2, 



and c are sets of four constants, 
J 



(7) 
(8) 
(9) 

(10) 



and the "corrector": 



1' J + 1 1, J - 1 3 1 i, j +1 ij 1. j - 1 



y 



j 3, 4, 5 . . . 



i 1. 2 N 



(11) 



Superscripts p and c indicate predicted and corrected values, respectively. 
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An RKG solution is obtained by successive applications of equations 
(3) through (6) for each j from to 3, together with four applications of 
equation (1). y are the values of the solutions at the mesh points determined 
byh. 

A Milne solution is obtained by using (10) followed by (1) to evaluate 
f^ • . 1 • Then (11) and another application of (1) determine y. . and f. . 
respectively. 

The maximum deviation (over i) between predicted and corrected values 
divided by 29 (see reference 1) is used to judge the size of the integration interval. 
The procedure for changing this interval will be described in detail below . 

Both of the above methods yield fourth order precision, i.e., have 
truncation errors of order h . 

THE PROGRAM 

Every subroutine, if it is to be complete and self-contained, requires of 
the user a set of instructions which enable the subroutine to perform its 
intended function and which supply to the subroutine the required return 
address(es) to the main program. Such a set of instructions is termed the 
"calling sequence . " 

This subroutine is designed in such a manner that the calling sequence, 
which is written as part of the user's program, has two distinct portions. 

The first portion is that in which the functions f. are to be programmed 
and stored in the appropriate locations . The data required for the computation 
of f j^ are available from given locations . 
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The second portion of the calling sequence is concerned with preparing 
the output and testing for continuation of the solution. As mentioned previously, 
both of these operations are entirely within the control of the user . The two 
parts of the calling sequence are separated by a TRSUB instruction. 

This separation was necessary because we do not wish to execute the 
output program every time control is transferred to the first part of the calling 
sequence. For example, when using the Milne method, y , y , y^ ', y 
are computed in that order. The second and fourth of the above quantities are 

computed in the first part of the calling sequence . The second part of the 

(c) 

calling sequence, however, is to be executed only after computation of y'^ ' . 

Three program parameters must be supplied by the user before execution 
of the program . 

1 . INI is entered into location 997 

2. h is entered into location 998 

3 . T is entered into location 999 

The initial conditions y , y „, y » . • .,y are entered into locations 500 

"^00 10 20 NO 

to 500+ INI . All input is in normalized floating point form. After the above 
data are stored, the calling sequence is as follows: 

loc . p : TRSUB +0 204 q 600 

loc. q : Start computation of f., storing them in 550 + i, i = 1, 2, . . .,N. 

loc.r: TRSUB +0 204 s 620 
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loc. s : Select output from Iocs. 500 to 500+lNI and 551 to 550+1 N| 
for punching. 
Punch output. 

Select from the same locations the criterion for terminating 
the solution. 
Test this criterion. 

If solution is to be continued, transfer to loc. 640. 
The y. . values for computing the functions f. between locations q and r are 
available from locations 500 to 500+ IN I . 

If the parameter N is stored as a positive number, the program will compute 
the first three points (not counting the initial point) with the RKG method. Then 
it will switch to the Milne method until the integration interval is to be 
decreased. If N is stored with a negative sign, the program will use the RKG 
method throughout. 

(If it should be desirable to compute m RKG points (m>3) before switching 
to the Milne method, this can be achieved by placing behind the program deck 
the instruction 

loc . 625: Set B +0 050 655 m + 1 .) 
The subroutine may be used any number of times within the same program . 
The only restriction here is that IN I is not changed. If iNJis to be changed, the 
program must first be reloaded. 

It might be pointed out here that the system of equations to be solved need 
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not be simultaneous . Several unrelated differential equations may be solved at 
the same time. 

DOUBLING THE INTERVAL 

When using the Milne method, the routine examines at each point the sum 

of the absolute values of the maximum deviations for the last three points . When 
this sum is exceeded by the quantity 29x10" ■'■/D, the interval is doubled. The 
value of D is largely a matter of experimentation. Obviously, there has to be 
some separation between the regions in which the interval is decreased or 
increased. The region in which the interval is left unchanged is given by 

29x10""^ ^ E < 29x10"'^. 
In fact, the ratio of these two bounds must be at least 30, since the error is 
multiplied by approximately that when the interval is doubled. It was found here 
in one specific instance that D = 200 was near a machine time optimum. 
However, D, which is located in 788, may be changed should the user so desire. 

After doubling of the interval has taken place, the routine is prevented from 
doubling again until 3 more RKG points, followed by one Milne point, have been 
computed . 

It has been found useful to indicate on the output any modification of h which 
may have taken place . Therefore the program will punch a card with a two- in 
the second word when the doubling routine is called in. 

DECREASING THE INTERVAL 

Before going into this process it is necessary to specify the two 
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precision criteria which were mentioned previously. 
If the entire program deck is used, the quantity 



E = max 
J i 



y(P) - y(^) 
ij ij 



will be compared to 29x10 ^ . If the last card is removed (Card 70 Deck 1), 
then the quantity 



E = max 
J i 



(P) 

1- ^ij 



y(9) 

ij 
-T 
will be compared to 29x10 . These are the fixed and floating point precision 

criteria, respectively. 

Three different routines for decreasing the interval are built into this 

program: 

1 . If the interval is to be decreased at point p in the middle of a run (and 

if h has remained unchanged at the previous point p ), the program will discard 



that point, back up to the previous point p , re-enter the RKG method and 
compute three points with the decreased interval . It will then switch again 
to the Milne method. An indicator card with an H will be punched. 

2. If, after having just decreased the interval, the program finds upon 
leaving the RKG procedure that it is to decrease the interval again, it will 
again back up to the original point p and proceed from there as before . The 
program will repeat this procedure as often as necessary. The results at 
point p are not repunched during this operation. 

3. If the point p happens to be the initial point of the solution, i. e., if 
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the initial h has been chosen too large, the program will, upon leaving the RKG 
procedure, punch a spacer card, repunch the initial point and start again with 
the reduced interval. It will repeat this procedure until the proper initial h 
has been found. 

If, in this last case, the programmed switch is set to STOP, a conditional 
stop will occur, displaying the value of h which caused it. This was done in 
order to give the user the opportunity to modify the starting procedure . If, 
when the stop occurs, the storage entry switch is set to plus and the start 
button is pressed, routine c will be executed. If the storage entry switch is 
set to minus, routine a will be executed. 

There is one more special situation which may arise . Suppose we are 
solving a single differential equation of large order N . Then it is conceivable 
that we may not wish all derivatives up to order N-1 to be examined for precision. 
Suppose we wish to examine only the first n (n < N) functions, i.e., 
y, y', y", . . ., y . This can be achieved by placing behind the subroutine 

a card loading into location 723 the instruction 

+0 109 B 717, 
where B =* 1000 - N + n. 

The fraction H is stored in location 789 and may be changed by the user, if 
he so desires . 

STORAGE 

Permanent. The program occupies locations 600-999. This breaks down 
as follows : The interpretive program (Deck 1), which corresponds to the 
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entire solution of the problem, occupies 180 locations. This interpretive 
program contains some forty-odd instructions which are functions of the 
program parameter I N I . Those instructions are compiled with two basic 
language programs (Deck 2) totalling 190 instructions . The basic language 
programs, which contain no loops, are executed only once. Numerical constants 
(Deck 3) take up 20 locations, and parameters 10 locations. 

Erasable. The buffer between the user's program and the subroutine, Iocs. 
500 to 500+1 N I and 551 to 550^-|NI , is erasable storage before and after 
execution of the subroutine . Location 550 cannot be used during execution of 
the subroutine. The remainder of these two bands is not used. 

If the RKG method only is used (N negative), the erasable storage requirement 
is locations 500-2(| N I +1) to 499 . 

If the Milne method is used (N positive), the erasable storage requirement 
is locations 500-13 |NI to 499. The appropriate block of erasable storage which is 
used to carry along the function values at previous points may not be used 
during execution of the subroutine, but is available before and after. For 
the sake of compactness the erasable storage was so arranged as to leave 
free as much of lower memory as possible . For all practical purposes N is 
unlimited in size if the RKG method is used. In the case of the Milne method, 
which requires much more storage, N = 30 leaves about 150 locations for the 
main program. The overall limit on N, dictated by the size of the buffer, is 49. 
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TIMING 

The routine was so designed that no operations are performed unnecessarily 

or repeated needlessly. N is tested only once, and the basic control operations 
are executed only once . It is felt that this routine is nearly as efficient as a 
custom-made program. If the RKG method is to be used, then no instructions 
pertaining to the Milne method are executed, with the exception of the one- 
time execution of a single SET C instruction at the beginning of the program. 

The execution time per point, excluding the time required for the 
computation of f, is about (6+3 I N I ) seconds for the RKG method and 
(2.5 + 1.5 I N| ) seconds for the Milne method. 

A detailed time study for a particular problem involving a single second 
order equation showed that the custom-made program required 14 seconds 
per point whereas the subroutine required 16 seconds per point. 

MISCELLANEOUS REMARKS 

The loop box and COUNT register are available anywhere for the user's 

own program. A word of caution is appropriate at this point. It is to be noted 
that the program in which the functions f are computed by the user is 
executed four times for each RKG point and twice for each Milne point. This 
is important if the coefficients of the equation are given in tabular form . It 
would then be necessary to count the number of executions before giving the 
loop instruction. 
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Two unconditional stops may occur during execution of the subroutine : 

1. A stop at location 635 indicates |Nl < 1 or INI > 100. 

2. A DIV-CHECK at location 717 indicates y^^^ = 

i, J +1 
for some i, if the floating point precision criterion is used. 

When tracing, the program is so arranged as to trace through each of the 
two sections of the user's program only once. Therefore it is advisable, when 
checking out this portion of the program, to enter the routine with some non- 
trivial initial conditions . 

The program may be translated with the standard Bell translation program. 

Deck numbers distinguish the three tjrpes of cards used. 

The error term E is stored in location 659 so that the skeptic may include 
J 

it in his output if he so desires". 

EXAMPLE 

The program for a specific sample problem is shown in Figure 1 . The 

equation to be solved is 

y^^^K y'y"y'"+ yy'= sinx, 
subject to the initial condition 

y(0) = 1, y'(0) = 0. y"(0) = 1, y'" (0) = 0. 
To reduce this equation to a system of first order equations, we define 
y' =s u 
u' = V 
v' = w. 
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Then w' = sin x - uvw - yu. 

The initial conditions become 

y(0) = 1, u(0) = 0, v(0) ^ 1, w(0) = 0. 
Four significant figures are required, and the initial h is taken to be . 1 . 

The output is to consist of the quantities x, y, u, v, w, and w' punched on 
one card. 

The solution is to be terminated when x exceeds unity or y becomes negative, 
whichever happens first . 

The entire program requires 17 instructions. 

FLOW DIAGRAMS 

Complete flow diagrams for this subroutine are shown in figures 2 and 
3 . A numeral or letter outside an instruction box denotes the location of the 
instruction represented by that box. 

Locations o^ and R contain the transfers to the main program , 
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SUMMARY TABLE 

In order to aid in the application of this program a summary table is given 

below : 

Item Purpose Location 



±N 
h 
T 

y, 







f. 
1 



1st address 
2nd address 
3rd address 

H 

D 

E. 
J 

+ 0050655 (m + 1) 
+ 0109B717 

Storage 
Storage 
Storage 
STOP TRACE 



Input 
Input 
Input 
Input - Output 
Input - Output 
Output 
TR to SUBR, initial . 
TR to SUBR, after computation of f . 
TR to SUBR, for continuation of 
solution 
Parameter 
Parameter 
Error term 
To compute m initial RKG points 
To test only up to y^""^^ 

B = 1000 - N + n 
Locations used by program 
Erasable storage for RKG 
Erasable storage for Milne 
Tracing control 



997 

998 

999 

500 
501 to 500+ INI 
551 to 550+ IN I 

600 

620 

640 
789 
788 
659 
625 

723 
600 to 999 

500-2(1 NI+ 1) to 499 
500-131 Nl to 499 
601, 621, and 641 
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FIGURE 1. PROGRAM FOR THE SOLUTION OF 

y^'VyV'y"'+yy'=sinx 
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9 
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+ 
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PUNCH 
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+ 


2 


500 


505 


000 


X-l 


113 


TRS 


+ 





201 
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114 


TEST X-l 


114 


MOV 


+ 


9 
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501 


000 


Y 


115 
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+ 





201 


640 
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+ 
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000 
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9 
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999 




500 




+ 









00 
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FIGURE 2. FLOW CHART FOR RUNGE-KUTTA-GILL METHOD AND 
STARTING PROCEDURE FOR MILNE METHOD 
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FIGURE 3. FLOW CHART FOR MILNE METHOD 
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DOUBLE TABLE LOOK-UP ON THE IBM 650 

Robert H. Goers s 
Dodco, Inc. 

INTRODUCTION 

One of the most difficult types of computation in a computing laboratory 
using digital equipment is double table look-up of location of data presented 
in a three-dimensional form. Specifically, Y is presented as a function of 
two variables X and Z. 



■—-— ^^^^^H^ 








^^^^*-^\r 


^Ns. Figure 1 


"^^i> ^^""^-^^ 


\^\. 


^-^ ^ ^^\ 


^Ns^^ N. 




^^ 



X 

For example, the presentation of drag of both aircraft and missiles is 

frequently in this form where drag is a function of velocity and altitude . 
One method of approach to this problem is to reduce the data to a function 
of one variable or a fit in either one or both variables . Unfortunately, 
this laborious procedure must be repeated for each set of data. 



Now associated with Electronic Associates, Inc, 
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This paper presents a general method of double table look-up of Y as 
a function of two variables X and Z . The procedure outlined has the advantage 
that tables of various sizes can be handled without changes in the code, 
subject to storage limitations of the equipment used. While the basic concept 
can be used on any medium or large-scale stored-program computer, this 
paper presents the storage limitations, code, and an example of a 5 x 5 
table for the IBM 650 . 

METHOD 

On the 650 the storage limitations are that the number of points in 
both X and Z be less than 48 and the number of values of Y be less than 
1800 minus the storage requirements of the problem using this subroutine. 

The subroutine will find the necessary values for linear interpolation 

for Y as a function of any X and Z in the range given by the table . The only 

restriction on the tabulated data is that Y must be given for the same points 

of X for all values of Z and vice versa. To interpolate linearly for Y as a 

function of X , Z , where X < X < X and Z -^ Z < Z , we need the following 
OZ iOIIAZb 

numbers: X^, X„, Z , Z„, Y,„ ^ ,, Y , Y , Y 

^ II A B (XjZ^) (Xj^z^) (XjZg) (XjjZg) 

i 
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To make the explanation of the presentation of the data and the code 
easily understood, we define X as the primary argument and Z as the 
secondary argument. The arguments and data are fed into the 650 as 
follows: The primary arguments are loaded in ascending order starting in 
storage location 0050; the secondary arguments are loaded in ascending 
order starting in location 0100; and the table of Y's are loaded in ascending 
order first with respect to the primary argument, then with respect to the 
secondary argument. Tables 1, 2 and 3 show the storage of primary and 
secondary arguments and the functions for a 5 x 5 table . 



TABLE 1 

Location of Primary Argument 
0050 X, 



0051 
0052 
0053 
0054 



X. 



TABLE 2 

Location of Secondary Argument 

0100 Z = a 

0101 Z = b 

0102 Z = c 

0103 Z= d 

0104 Z = e 



TABLE 3 
Storage of Function Y(XZ) 



x= 1-5 



z = a-e 



Location 


Contents 


Location 


Contents 


Location 


Contents 


0200 


Y(la) 


0210 


Y(lc) 


0220 


Y(le) 


0201 


Y(2a) 


0211 


Y(2c) 


0221 


Y(2e) 


0202 


Y(3a) 


0212 


Y(3c) 


0222 


Y(3e) 


0203 


Y(4a) 


0213 


Y(4c) 


0223 


Y(4e) 


0204 


Y(5a) 


0214 


Y(5c) 


0224 


Y(5e) 


0205 


Y(lb) 


0215 


Y(ld) 






0206 


Y(2b) 


0216 


Y(2d) 






0207 


Y(3b) 


0217 


Y(3d) 






0208 


Y(4b) 


0218 


Y(4d) 






0209 


Y(5b) 


0219 


Y(5d) 
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Inspection of Table 3 indicates that it could be considered a table of 
tables; that is, tables of Y as a function of X are arranged as in 
ascending order of Z. Associated with both arguments is a table of 
secondary modifiers . The secondary modifiers serve as an index to 
which table of Y as a function of the primary argument should be used. 
Specifically, adding the data address of the secondary modifier to the 
location of Xj yields the location of Y. . . 

TABLE 4 
Storage of Secondary Modifiers 



Location 


Contents 




0150 


00 0150 0002 


z = a 


0151 


00 0155 0002 


z = b 


0152 


00 0160 0002 


z = c 


0153 


00 0165 0002 


z = d 


0154 


00 0170 0002 


z = e 



Associated with the primary argument is a number called the jump 
modifier . Adding the data portion of the jump modifier to the data address 
of Yjj^ yields the data address of Y . It is important to note that for a 
table of any size the only changes that need to be made in the routine other 
than the data and arguments are the secondary and jump modifiers . 

The code itself consists of 31 instructions, which logically break 

down into eig^t steps . In each of these steps a Store Distributor command 

with a variable data address is generated in the lower accumulator . This 

operation is then executed, and the contents of the distributor are stored 

in a desired storage position . The arguments X and Z are assumed 

u ^ 

stored in locations 0030 and 0031 . The following operations are performed: 

1 . Find the location of Z by 650 table look-up; store Z in 0035- 

B 6 

2. Find the location of Z by subtracting 1 from the data address of 

Z^; store Z in 0034. 
° A 
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3 . Find the location of the secondary modifier by adding 50 to the 
address of Z.; store the secondary modifier in 0046. 

4. Find the location of X by 650 table look-up; store X_ in 0033. 

5 . Find the location of X by subtracting 1 from the location of X ; 

store X in 0032. 
I 

6 . Find the location of Y by adding the secondary modifier to the 
location of X ; store Y in 0036. 

7 . Find the location of Y,^ . by adding 1 to the location of Y ; store 

IIA ' s lA 

Y„ , in 0037 . 
IIA 

8 . Find the location of Y by adding the jump modifier to the location 

of Y„ . ; store Y^ , in 0038 . 
IIA lA 

9 . Find the location of Y by adding 1 to the location of Y ; store 

ilB IB 

Y^gin0039. 
The complete code and a list of storage positions follows: 



TABLE 5 
Storage of Input and Output 

Location Contents 



0030 Xq 

0031 Z^ 

0032 Xj 

0033 Xji 

0034 Z^ 

0035 Zb 

0036 Yj^ 

0037 Yjj, 

0038 Yjg 

0039 Y„3 

Jump Modifier 

0040 00 0004 0002 
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TABLE 6 

Other Storage 
Location Contents Function 

0041 00 0000 9998 Address Modifier 

0042 00 0050 0002 Address Modifier 

0043 00 0001 0002 Address Modifier 

0044 69 0000 0004 Load Distributor Command 

0045 69 0000 0012 Load Distributor 

0046 Temporary Storage 

This system was used by the author as a subroutine in a problem 
coded with the floating decimal interpretive system described by 
V. M. Wolontis in IBM Applied Science Division Technical Newsletter 
No. 11 . For that reason the data address of the last command is 1095, a 
location for a return to the Wolontis system . In this connection it should 
be pointed out that when using any floating coded decimal system, or when 
using a floating decimal arithmetic imit, care must be taken that the 
argument is only one order of magnitude . This can be accomplished by 
using a modified argument equal to the argument increased by the quantity 
1 X 10 oC, where oCis such that all modified arguments have the same 
exponent . 

If there are n values of X and m values of Z, the secondary modifiers 
and jump modifier can be expressed in terms of n and m as follows: 

Command Data Instruction 

Jump modifier 00 

Secondary modifiers 00 

00 

00 

00 



00 0150 + (m-l)n 0002 



n - 1 


0002 


0150 


0002 


0150 + n 


0002 


0150 + 2n 


0002 


0150 + 3n 


0002 
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With slight modifications, this system can be expanded for higher 

order interpolation or to handle two related functions of the same two 

arguments . 
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IBM IBM TYPE 650 PROGRAM SHEET 

Ttn-lc Mark 

PRORI FM- DOUBLE TABLE LOOK-UP WRITTEN BY: R. H. GOERSS 



FORM NO. 2Z-6I8I-I 
PRtNTEO IN U.S.A. 



INSTR 
NO. 


LOCATION OF 
INSTRUCTION 


OPERATION 


ADDRESS 


REMARKS 


ABBRV. 


CODE 


DATA 


INSTRUCTION 




0001 


RAL 


65 


0044 


0002 


Reset lower accumulator basic load 














distributor command. 




0002 


LD 


69 


0031 


0003 


Load Z^ (secondary argument) into the 














distributor . 




0003 


TLU 


84 


0100 


8002 


Find the location of Zg. 




8002 


LD 


69 


L (Zg) 


0004 


Load Z_ into distributor. 




0004 


STD 


24 


0035 


0005 


Store \ 0035. 




0005 


SL 


16 


0041 


8002 


Find the location of Z^ by subtracting 














1 from the location of Zr, add 2 to 














the instruction address. 




8002 


LD 


69 


L (Za) 


0006 


Load Z^ into the distributor. 




0006 


STD 


24 


0034 


0007 


Store Z^ in 0034. 




0007 


AL 


15 


0042 


8002 


Find the location of the secondary 














modifier and increase instruction 














address by 2. 




8002 


LD 


•69 


L (SM) 


0008 


Load secondary modifier into 














distributor. 




0008 


STD 


24 


0046 


0009 


Store secondary modifier in temporary 














storage. 




0009 


RAL 


65 


0045 


0010 


Reset lower accumulator load basic 














load distributor command. 




0010 


LD 


69 


0030 


0011 


Load primary argument Xq into the 














distributor. 




0011 


TLU 


84 


0050 


8002 


Find the location of Xj^. 




ftOO? 


T.n 


69 


T. rx,,"* 


0012 


Load Xjj into distributor. 




nm? 


STD 


7U 


00 ■^^ 


0013 


Store Xjj in 0033. 




00 Tl 


SL 


16 


0041 


8002 


Find location of Xj by subtracting 1 














from the location of X-j-^. 




fi002 


T.n 


fiq 


T. rxj.^ 


0014 


Load Xj into distributor. 




0014 


STD 


24 


0032 


0015 


Store Xj in 0032. 




001 s 


AL 


IS 


0046 


8002 


Find the location of Yj^ by adding 














the secondary modifier^to the location 














of Xj." 




8002 


LD 


69 


L (Yj^) 


0016 


Load^Y-j-j^ into the distributor. 




0016 


STD 


24 


0036^^ 


0017 


Store YiA in 0036. 




0017 


AL 


15 


0043 


8002 


Find the location of Yjjn by adding 














1 to the location of Yja* 




8002 


LD 


69 


L ^^T\h\ 


0018 


Load Ytt. in distributor. 




0018 


STD 


24 


0037 


0019 


Store Yjj^. 




001 q 


AT. 


t; 


0040 


8002 


Find the'"location of Y,„ by adding 














■jump modifier to the location of Yjj^. 




ft009 


T.n 


fig 


T. (y \ 


0020 


Load Y-j-jj into distributor. 




0090 


STn 


24 


00.38^^ 


0021 


Store Y^g in 0038. 




0021 


AL 


15 


0043 


8002 


Find location of Y„_. 




R007 


T.n 


fiQ 


T. (Ytth) 


0022 


Load Yjjg in distributor. 




0077 


STD 


24 


0039 


1095 


Store tTTTt- 
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